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ABSTRACT 

Collapse of the rotating magnetized molecular cloud core is studied with the 
axisymmetric magnetohydrodynamical (MHD) simulations. Due to the change of 
the equation of state of the interstellar gas, the molecular cloud cores experience 
several different phases as collapse proceeds. In the isothermal run-away collapse 
(n ^ 10^°H2cm~^), a pseudo-disk is formed and it continues to contract till the 
opaque core is formed at the center. In this disk, a number of MHD fast and 
slow shock pairs appear running parallelly to the disk. After the equation of state 
becomes hard, an adiabatic core is formed, which is separated from the isothermal 
contracting pseudo-disk by the accretion shock front facing radially outwards. 
By the effect of the magnetic tension, the angular momentum is transferred 
from the disk mid-plane to the surface. The gas with excess angular momentum 
near the surface is finally ejected, which explains the molecular bipolar outflow. 
Two types of outflows are observed. When the poloidal magnetic field is strong 
(magnetic energy is comparable to the thermal one), a U-shaped outflow is formed 
in which fast moving gas is confined to the wall whose shape looks like a capital 
letter U. The other is the turbulent outflow in which magnetic field lines and 
velocity fields are randomly oriented. In this case, turbulent gas moves out 
almost perpendicularly from the disk. The continuous mass accretion leads to 
the quasistatic contraction of the first core. A second collapse due to dissociation 
of H2 in the first core follows. Finally another quasistatic core is again formed 
by atomic hydrogen (the second core). It is found that another outflow is ejected 
around the second atomic core, which seems to correspond to the optical jets or 
the fast neutral winds. 

Subject headings: ISM: clouds — ISM: magnetic fields — stars: formation — 
ISM: jets and outflows 
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1. Introduction 

Star formation has been a long-standing target in astrophysics. The infrared protostar 
distribution revealed that the molecular cloud cores, which coincide with relatively high- 
density part (n ~ lO^cm"^) of the molecular clouds, are the sites of star formation. The 
observed molecular cloud cores are divided into two categories: those observed associated 
with and without protostars. The molecular cloud cores without protostars are called starless 
cores or prestellar cores and are considered younger than the cores associated with protostars 
(protostellar cores). From a theoretical point of view, clouds or cloud cores experience the 
isothermal run-away collapse first and then accretion on to the stellar core develops (Larson 
1969). Prestellar cores which indicates inflow motions (e.g. in L1544 rotation and infall 
velocities ~ 0.1km s~^ are observed by Ohashi et al. (2000)) do indicate that they are in 
the dynamically contracting phase, in other words, in the run-away collapse (Ciolek & Basu 
2000). After the epoch when the dust thermal emissions are trapped in the central part 
of the cloud {ric ~ lO^^cm"^), an adiabatic core is formed and isothermal gas continues to 
accrete on to the core. The molecular cloud core in this phase is observed as a protostellar 
core. It is shown that the dynamical evolution of the cloud core is characterized by the 
sequence from the prestellar cores to protostellar cores. 

Dynamical collapse of the magnetized clouds are studied by many authors (Scott & 
Black 1980; PhiUips 1986a,b; Dorfi 1982, 1989; Bentz 1984; Mouschovias & Morton 1991, 
1992; Fiedler & Mouschovias 1992, 1993; Basu & Mouschovias 1994; Tomisaka 1995, 1996; 
Nakamura et al. 1999). Rotating clouds collapse has been attacked seriously with numerical 
simulations (Bodenheimer, Tohline, & Black 1980; Norman, Wilson, & Barton 1980; Wood 
1982; Narita, Hayashi, & Miyama 1984; Truelove et al. 1997, 1998; Tsuribe & Inutsuka 
1999a,b). However, a restricted number of articles are published regarding the dynamical 
contraction of the cloud with both rotation and magnetic fields (Dorfi 1982, 1989; Basu 
& Mouschovias 1994, 1995); for quasistatic evolution see Tomisaka, Ikeuchi, & Nakamura 
(1990). These researches are confined to the relatively early prestellar stage. 

Is it sufficient to consider the effects of the rotational motion and the magnetic fields 
separately? In the dynamical contraction phase, it is shown that the molecular outflow is 
driven by the cooperative effect of the magnetic fields and rotation motion (Tomisaka 1998). 
The toroidal magnetic fields are generated from the poloidal ones by the effect of rotation 
motions. The magnetic torque works only when the poloidal and toroidal magnetic fields 
coexist. The magnetic torque and thus magnetic angular momentum transfer along the 
magnetic field line is important to eject the outflow. Since the outflow brings the excess 
angular momentum, the angular momentum that remain in the adiabatic core and thus 
new-born star is reduced by a factor from 10"^ to 10~^ from that of the parent molecular 
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cloud core (Tomisaka 2000). The outflow have not observed either in the magnetized but 
non-rotating cloud (Scott & Black 1980; Tomisaka 1996) or the rotating but non-magnctized 
cloud (Norman, Wilson, & Barton 1980). Therefore, rotation and magnetic fields are both 
essential to the evolution of molecular cloud cores. In the present paper, we present the 
dynamical contraction of the magnetized and rotating cloud. 

The cooperative effect of magnetic fields and the rotation motions becomes important 
after the adiabatic core is formed near the center of the cloud core (Tomisaka 1998). There- 
fore, the evolution throughout from the prestellar to protostellar core should be studied. 

Plan of this paper is as follows: in §2 model and numerical method are described. As 
the initial condition, we choose a slowly rotating cloud with purely poloidal magnetic fields 
(no toroidal magnetic fields). And we follow the evolution using magnetohydrodynamical 
simulations. Section 3 is devoted to the numerical results. In this section we compare clouds 
with strong magnetic fields and those with weak magnetic fields. This shows that completely 
different two types of outflows are ejected in respective clouds. Another comparison is made 
between fast rotators and slow rotators. In §4, we discuss the evolution till the second core, 
which becomes actually a new-born star, is formed. It is found that another outflow is found 
around the second core, which seems to correspond to the optical jets or high speed neutral 
winds. We also discuss whether the mass inflow/outflow rate and the momentum outflow 
rate observed in molecular bipolar outflows are explained or not. 



2. Model and Numerical Method 

To study the dynamical contraction, we consider an isothermal cylindrical cloud in 
hydrostatic balance with inflnite length as the initial state. In terms of the gravitational 
potential ijj and isothermal sound speed Cg, the radial distribTitions of density p, rotation 
speed v^, and magnetic field density are calculated using the equation of hydrostatic 
balance as 

vl cldp 1 dBl 



r dr p dr Sirp dr 
and the Poisson equation for the self-gravity as 



- 0, (1) 



where G represents the gravitational constant. This equation has a solution as follows 
(Stodolkiewicz 1963) 

Po{r) = Pco (^1 + , (3) 
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/ ^2 X -1/2 

and 

B.cW = B.„(^)"^iS.„(l + ^) (5) 



with using the scale-height at the center, 

H! = ?;."^"^:"r . (6) 



where Pco, ^co, ^-nd B^o represent respectively the density, the angular rotation speed, and 
the magnetic flux density at the center of the cylindrical cloud (r = 0). This is the same 
hydrostatic configuration studied by Matsumoto, Hanawa, & Nakamura (1997). However, we 
assume no initial toroidal magnetic field -B,^ = in contrast to them. The density distribution 
extends till the radius where the thermal pressure c^ps becomes equal to the external pressure 
Pext- The solution contains three non-dimensional parameters characterizing the distribution 
after adopting natural normahzation such as for distance r' = r/H = r /[cs/{AnGpsy^'^], 
time t' = t/rs = t/[l/(47rGp5)^/^], and density p' = p/ps = c^gp{r) / Pe^^. We summarized 
the conversion factors from non-dimensional to physical quantities in Table 1. The first 
parameter characterizing the initial state is related to the magnetic to thermal pressure 
ratio as 

a = B,{rf/A7^p{T)cl = BlJAi^p^^cl (7) 
and the second one is related to the angular rotation speed as 

il' = iU/{ATiGpsY'\ (8) 

Finally the surface to center density ratio 

F = Peo/Ps- (9) 
The scale-height at the center He is expressed using these non-dimensional parameters as 

_ l + a/2 

Prom this equation, it is shown that a hydrostatic balance is achieved only when F > 
2Q'^, i.e., AnGpco > Sil^o- Matsumoto, Hanawa, & Nakamura (1994)'s /3 parameter, which 
represents the ratio of the centrifugal force to the thermal pressure force, is expressed using 
our parameters as 

^ ^ p{r)v<p{rf/r 
c^dp{r)/dr 
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2Q 



,2 1 + «/2 
F - 2f]'2 ■ 



(11) 



To initiate the cloud collapse, we added a small density perturbation. The wavelength 
of the perturbation is taken equal to that of the most unstable Jeans mode. The linear 
instability of the rotating, magnetized, isothermal cylindrical cloud is studied by Matsumoto, 
Hanawa, & Nakamura (1994). They gave a fitting formula for the most unstable wavelength 
as follows: 



A. 



(47rG'p,)V2 



-1 



27r(l + a/2 + /5)'/' 



0.72 



;i + «/2 + /5)^/^ - 0.6 
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using equations (7), (9), and (11). Therefore, we take this most unstable mode and the 
initial density is assumed equal to 



p{z,r) = po{r) 



1 + Acos 



27rz 



(13) 



for — Amax/2 < z < Ainax/2. The amplitude of the perturbation A is taken equal to 0.1. 
Hereafter, we omit the primes which indicate the non-dimensional quantities, unless the 
quantities are confused with dimensional ones. 

Prom calculations of one-dimensional, non-rotating, non-magnetized cloud, the effective 
equation of state of the gas forming a star is summarized as follows (Tohline 1982). Gas in 
the molecular cloud core with p ~ 10^H2 cm^'^ obeys the isothermal equation of state with 
the temperature of Tq ~ lOK. However, the molecular gas becomes optically thick against 
the thermal emission from the dust, which cools the molecular core, after the contraction 
proceeds and the central density reaches Pc = Pa ~ 10^°H2 cm~^ (Masunaga & Inutsuka 
1999). The gas becomes adiabatic beyond the density and the first core forms (Larson 
1969). Hydrogen molecules (H2) dissociate into atomic hydrogen (H) when its temperature 
reaches Tdis — lO'^K. Typical density at which the temperature reaches Tdis is equal to 
Pb = Pa (Tdis/To) 1/(^-1) ~ lO^^Hs cm-3(Tdis/10^K)5/2(To/10K)-5/2. The dissociation ends at 
the density pc — 10^^H2 cm^'^. Finally forms an atomic hydrogen core, which is call the 
second core by Larson (1969). To include these changes, we adopt the multiple polytropes 



P 



7/5 



P < PA, 

Pa < p < Pb, 

Pb < P < PC, 



(14) 
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with PA = lO^^Hscm-^ = lO^^Hscm-^ and pc = lO^'^Hacm-^ 



The basic equations to be solved are the magnetohydrodynamical equations and the 
Poisson equation for the gravitational potential. In cylindrical coordinates {z, r, 0) with 
d/dcf) — Q, the equations are expressed as follows: 



dp d . , 19, , 



0, 



(15) 



dpvz d , .Id. . 



,dp 



dilj 
dz 



"^Wz-P-^- 



1 

47r 



dz 



-Brh — 



dB,. dB, 
dz dr 



Br 



(16) 



dpVr 

"df 



+ 



d , ,19, , 

^[pVrV^) + - — [rpVrVr) 

dz r dr 



'dp 



dip 
dr 



1 

47r 



~|:(rs,)s,+ 



dBr dB, 



dz dr 



B. 



dprv^ d N , 1 ^ / N ^ 

+ g-^{prv,v,) + -g-^irprv,vr) = - 



'-^irB,)B^ + ^B. 
r dr dz 



(17) 
(18) 



dB^ 
dt 
dBr 

~dr 



ld_ 
r dr 



[r{v^Br - VrB^)], 



d 

- — {v^Br - VrB^), 

dz 



dB^ 
dt 



v.B^ 



d_ 

dr 



VrB^ - V^Br), 



d'^ip 1 d 

— - H 

J. Qj. 



dip 
dr 



AnGp, 



(19) 
(20) 
(21) 
(22) 



where the variables have their ordinary meanings. Equation (15) is the continuity equation; 
equations (16), (17) and (18) arc the equations of motion. The induction equations for the 
poloidal magnetic fields arc equations (19) and (20) and for the toroidal magnetic field is 
equations (21). The last equation (22) is the Poisson equation. 

The MHD equations were solved using van Leer (1977)'s monotonic interpolation and the 
constrained transport method by Evans & Hawley (1988) with the method of characteristics 
(MOC) modified by Stone & Norman (1992). To ensure the specific angular momentum prv<^ 
and the toroidal magnetic fields 5^ are convected similar to the density, advections of such 
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values are calculated using the consistent advection which is first pointed out by Norman, 
Wilson, & Barton (1980). We rewrote equations (18) and (21) as 





d 


dt 


dz 


dprv^ 


d 


dt 


dz 



_ d_ 
dr 



(23a) 
(23b) 



To evaluate the right-hand-side of equations (18) and (21), we chose the same numerical 
mass flux pv used in equations (15) and (16) and multiplied B^/p and rv^p to the mass 
flux to get numerical flux for the angular momentum density and the toroidal magnetic 
fields. To solve the Poisson equation (22), we adopted MILUCGS [modified incomplete 
LU decomposition preconditioned conjugate gradient squared method: Meijerink & van der 
Vorst (1977); Gustafsson (1978)]. 

The gravitational contraction proceeds in a significantly non-homologous way. Thus, 
to see the late phase evolution, high spatial resolution is needed especially near the center. 
Similarly to the previous papers (Tomisaka 1996, 1998, 2000), we adopt the nested grid 
method (Berger & Oliger 1984; Berger & Colella 1989). In this method, a number of levels 
of grids with different spacings are prepared; finer grids cover the central high-density portion 
and the coarser ones cover the cloud as a whole. The grids are named as LO (the coarsest), 
LI, L2, ... and the grid spacing of Ln is chosen equal to a half of that of Ln — 1, that is, the 
grid spacings of Ln is equal to Azn — A.zo2~'" and Ar„ = Aro2~". The LO grid covers the 
region -Amax/S < z < and < r < Amax- 

In the nested grid method, the boundary condition is applied to the outer boundary of 
the coarsest grid LO, on which we adopted the fixed boundary condition at r = Amax and the 
periodic boundary condition at 2; ± Amax/S. The code has been tested comparing the results 
obtained with and without the nested grid technique [for detail see Tomisaka (1996)]. The 
the most unstable growth rates of perturbations are slightly different (~ 8%): pc = 10^ was 
attained at t = 1.25 for without nested grid (but using relatively large number of grid points 
400x400), but at t = 1.351 in the calculation with the nested grid method. This is mainly 
due to the fact the eigen-mode of the most unstable mode is as wide as the whole size of the 
numerical box and this mode is captured better by a simple calculation with a large number 
of zones (400x400) rather than the LO grid (64 x 64) of the present scheme. However, after 
the perturbation has grown to be non-linear, the growth rate agrees well. 
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3. Results 

3.1. Dynamical Contraction 

In Model A, we calculated the evolution with a — 1, Qq — 5, F — 100, and Ps — 
10^H2cm~^. We summarized the adopted parameters in Table 2. Similar to the non- 
rotational magnetized cloud [see Figs. 2b and 2c of Tomisaka (1996)], the cylindrical cloud 
fragments into prolate spheroidal shape whose wavelength is equal to Amax- This prolate 
spheroidal shape coincides with the structure expected from the linear stability analysis by 
Matsumoto, Hanawa, & Nakamura (1994). Next, this density enhanced region begins to 
contract along the major axis of the cylindrical cloud, since the magnetic fiekls are assumed 
to run parallelly to the major axis. Finally it forms a contracting disk (pseudo-disk) perpen- 
dicular to the magnetic field lines [Fig. 2d of Tomisaka (1996)]. The snapshot at this stage 
is shown in Figure 1. Using the conversion factor shown in Tablel, the epoch t — 0.6066Tff 
corresponds to 1.06Myr(ps/10^H2 cm~^)~^/^ from the beginning of calculation. Respective 
panels of Figure 1 have different spatial coverage. Figure la, which shows LI, captures global 
structure of contracting disk extending horizontally which is perpendicular to the magnetic 
field lines. The spatial resolution of LI is so limited that there seems no internal structures 
in the contracting disk. However, L5 which has 16 times finer resolution than LI shows shock 
fronts facing outer directions extend parallel to the r-direction. Fronts near z ~ ±0.02i/ are 
the fast-mode MHD shock front, because the magnetic fields bend toward the front passing 
the shock front. We can see another density jump near z ~ O.OIH (hereafter we will omit 
the sign ± and mention the upper half part of the figure since the structure is symmetric) . 

The shock fronts parallel to the disk are known in non-magnetized rotating isothermal 
clouds (Norman, Wilson, & Barton 1980; Matsumoto, Hanawa, & Nakamura 1997). This is 
not due to the rotation; multiple shock fronts are also found in the contracting magnetized 
cloud without rotation (Nakamura et al. 1999). However, situation becomes a bit comphcated 
in this cloud. In the outer region z ^ 0.02H, the magnetic field lines run almost vertically 
(Bz ^ Br and B^). Passing the MHD fast shock, in the intermediate region {O.OIH ^ 
z <^ 0.02H), the radial and toroidal components are amplified and the density increases 
compared to the outer region. Finally after passing another front near z ~ O.OIH, the 
toroidal component B^ decreases. That is, since the magnetic field lines defiect departing 
from the front at the second front, this is the slow MHD shock. The density range captured 
in this panel is from lO^^s = 10^-3(p,/100H2 cm'^^) to lO^^s = 109-3(p,/100H2 cm"^). 
It should be noticed that these shocks occurs in the isothermal gas. This phase is called 
"run-away collapse," in which the central density (pc) increases greatly in a finite time-scale. 

Figure Ic shows the structure captured by LIO. Almost all the gas in this figure is 
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isothermal. However, a central small part of the contracting disk r ^ 6 x 10~ if and 

l^l ^ 1 X 10~^H enters the adiabatic regime. At this stage, we can see another density jump 
is forming just outside the adiabatic part of the disk. It seems to grow into an accretion shock 
front, since it is known that an accretion shock forms outside the core when the adiabatic 
core develops (Larson 1969). This is easily understood as follows: the adiabatic gas with 
specific heat ratio 7 > 4/3 has a hydrostatic equilibrium irrespective of its mass. The scale- 
height in the z-direction of the core becomes larger than that of the disk. The adiabatic part 
of the contracting disk forms a spherical static core. 

Figure 2 shows the cross-cut view along two axes (panel a: along the disk mid-plane 
z — and panel b: along the z-a.xis r = 0). The lines with 0.6066 represent the stage 
shown in Figure 1. At this stage {t — 0.6066rff), inflowing gas is almost isothermal (p < pa)- 
In Figure 2a, the radial distributions of the density, the magnetic flux density, and the ra- 
dial and toroidal components of velocity are shown. We can see that density and mag- 
netic flux density distributions are approximately expressed by power-laws as p oc r"^ 
and Bz oc except for the central part. At t = 0.6067rfj, the density in the core ex- 
ceeds lOVs = 10^^(ps/100H2cm-^)H2cm-3. And at t = 0.6068rff, it reaches lO^V^ = 
10^^(ps/100H2 cm~^)H2 cm~^. At this stage, a radially outward-facing shock front is seen 
even inside the disk; infall motion is abruptly decelerated and the density and magnetic flux 
density are amphfied. This shows that a compact core is formed inside the accretion shock 
front. The central density increases with time and inside r ^ 1.8 x 10~^i? adiabatic gas (p > 
Pa) distributes. The size of the core is equal to r ~ 1.9xl0^^(cs/190ms~-^)(p5l0^H2 cm"^)"^/^. 
This reduces with time since the mass of the core increases by the effect of continuous ac- 
cretion. 

Before the shock front is formed t < 0.6067Tfr, the radial inflow velocity takes the 
maximum about ~ 2.5cs near r ~ 7 x 10~^H. For the Larson- Penston self-similar solution 
for the spherically symmetric dynamical collapse (Larson 1969; Penston 1969), this maximum 
inflow speed is expected equal to ~ 3.28cs. On the other hand, it equals to ~ 1.736cs for 
non-rotating isothermal disk (Saigo & Hanawa 1998). Therefore, it is shown that the actual 
inflow speed ranges between the spherically symmetric self-similar solution and the axially 
symmetric thin disk solution. After the shock front is formed around the core, the inflow 
velocity takes the maximum just outside the shock front and the maximum speed increases 
with time. Inflow motion is accelerated toward the shock front. Similar acceleration is also 
seen in the toroidal velocity, t),^. Before the core formation the toroidal speed v^j) takes the 
maximum near r ~ 2.5 x 10~^H. However, increases toward the accretion shock after 
core formation. At t — 0.6068Tff it reaches ~ Scg (see also Fig.l of Tomisaka (1998)). 

Structure seen in the cross-cut along the 2;-axis is more complicated (Fig. 2b). Two shock 
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fronts mentioned earlier (Fig. lb) correspond to the jumps near l^j ~ 0.02H and ~ O.OOSif^ 

At t = 0.6066rfj, the density and the inflowing velocity distributions have no discontinuities 
besides these two shock fronts. However, at t = 0.6067rfj, the inflowing velocity distribution 
begins to indicate a clear discontinuity near \z\ ^ 1.5 x lO^^H. This is a newly formed shock 
front and propagate spatially. Comparing two curves of t = 0.6067rff and t = 0.6068rfj, it is 
shown that this shock front breaks into two fronts and the inner one {\z\ ~ 1.5 x 10^"^ if) is 
standing still, while another outer one {\z\ ~ 5 x lO'^^H) is propagating outwardly. These 
two shock fronts are outwardly facing. Thus the inwardly propagation of the inner fronts is 
due to the infalling gas motion. 



3.2. Outflow 

Figure 3 illustrates the structure at t = 0.6069rflf. Although the gas is inflowing both 
inside and outside of the disk ai t — 0.6067rg (Fig. 2), at this stage t — 0.6069rff (r — 3.2 x 
10~^Tff)^ prominent outflow is formed. It is shown that the flow pattern is completely changed 
in At ~ 2 X lO^^Tff ~ 400yr. Outflow sweeps a sphere with radius of r ^ 1.2 x 10~^H (Fig.3a). 
Figure 3b indicates that the gas near the disk surface flows inwardly for r ;^ 2 x 10~'^H. 
However, the direction of the flow is changed upwardly near r ~ 2 x 10-^ H. Finally this 
gas is ejected. While the gas near the mid-plane of the disk {\z\ ^ 1 x lO^^H) continues to 
contract. This is reasonable because the total amount of angular momentum in one magnetic 
flux tube must be conserved in the axisymmetric ideal MHD simulation; for the outflow gas 
to get angular momentum, a part of the gas in the same magnetic flux tube has to lose its 
angular momentum and falls further. In the acceleration process of the gas, the angular 
momentum is transferred from the gas near the mid-plane to the gas near the surface of the 
disk. Considering the angular rotation speed, the angular momentum is transferred from 
the fast-rotating mid-plane to the slowly rotating surface gas. 

From Figures Ic and 3a (both show the structure captured by LIO), we can see the 
magnetic fleld lines run completely differently comparing before (pc < Pa' Fig-lc) and after 
(pc > Pa' Fig.3a) the adiabatic core formation. That is, in the isothermal runaway collapse 

phase (Fig.lc) the magnetic fields lines run vertically, in other words, perpendicularly to the 
disk. In contrast, after the adiabatic core is formed, the disk continues to contract and drags 



^Thc shape of the inner slow MHD shock is concave. On the z-axis it is found near z ~ 0.005H, while 
departing from the z-axis (r ^ 0.02_ff) it is fomid near z ~ O.OIH. 

"^t represents the time from the beginning of calculation but t represents the time after the core formation. 
We assumed that the core consists of the gas with density p > pA- 
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the magnetic field lines inwardly. Thus the angle between the magnetic field lines and the 
disk decreases. 

Figure 3b is a close-up view whose spatial resolution is 4-times finer than that of Figure 
3a. This panel shows us that the angle between the fiow and the disk is about ~ 45deg. 
The reason why the outfiow begins after the core formation is related to the angle between 
the magnetic field lines and the disk, ^mag- Blandford & Peyne (1982) have pointed out 
that for a cold gas rotating with the Keplerian speed to get angular momentum from the 
Keplerian disk via infinitely strong magnetic fields, ^^mag must be smaller than a critical 
value Oct = 60 deg. This is understood as follows: Consider the gas on one magnetic flux 
tube. When the magnetic flux tube rising steeply from the disk as 6'mag > the gas has 
to climb the effective potential well even if it rotates with the same angular speed of the 
Keplerian disk. On the other hand, when ^mag < ^cr, gas can escape from the gravitational 
well by getting angular momentum from the disk, if the gas has the same angular speed of 
the Keplerian disk. Although the exact value of 9cr depends on the disk rotation speed and 
the disk-to-central star mass ratio, small angle is preferable to acceleration. It should be 
noted that this configuration is achieved only after the core formation (Tomisaka 1998). 

Figure 3c shows the ratio of the toroidal magnetic component to poloidal one by contour 
lines, which is overlaid on the magnetic field lines. Since the toroidal-to-poloidal ratio is as 
small as ~ 0.6 in the disk, the disk is poloidal-dominated. However, in the region where the 
gas fiows outwardly the toroidal component grows and the toroidal-to-poloidal ratio reaches 
;^ 5 — 8. The outflow gas is toroidal-field-dominated. The coincidence of acceleration region 
with the toroidal-dominant region seems to indicate that the toroidal fields play an important 
role to accelerate the gas. The toroidal component of the Lorentz force. 



works below this toroidal-dominant region, that is, 2; ^ 5 x 10'^ H. This toroidal component 
accelerates the toroidal velocity and resultant toroidal motion amplifies the toroidal 
component of the magnetic fields. Outfiow speed exceeds the sound speed and the fastest 
speed reaches Vout — 7.5cs at this time. It increases with time. 



UzBr - jrBz) 




(24) 



3.3. Effect of the Hardness of the Polytropic Gas 



Although the outflow seems to continue, the further evolution is hard to study, because 
the time-scale (the free-fall time-scale at the central core) becomes shorter and shorter. 



- 12 - 



Therefore, we study model AH with a constant polytropic index larger than that of model 
A. In models AHl (Fig. 4a) and AH2 (Fig. 4b), the polytropic indices are chosen F = 2 and 
F = 5/3, respectively, for p > pA- (Models whose name have "H" have simple polytropic 
relation with F = 2orF = 5/3forp> pA-) Due to the hard polytropic index, the size of the 
adiabatic core, whose surface is determined by the jump in Vr, becomes large; for example 
at t = 0.6069rfr the size is equal to Tc — 4 x 10~*H for model AHl and Tc — 3 x 10~*H for 
model AH2 while it is Tc ~ 1.5 x 10~^H for model A. Similar to model A, just outside the 
core, outflow begins to be accelerated. The region which the outflow sweeps expands and 
the surface which separates inflow and outflow forms another MHD shock front. And the 
expansion of the front is very similar to that of model A (the front reaches z ~ 1 x 10~^if 
at this time which is similar to model A). 

To see the similarity in more detail, we calculated the mass of the core 



for models A, AHl, and AH2. These are equal to O.U07cy {ATiGf^'^ pl^^ , 0.1504c^/(47rG')3/V^^^ 
and 0.14504/ {4nGf/^pl^^ at the time t = 0.6069rff (r = 3.2x IQ-Vg) for models A, AHl and 
AH2, respectively At later epoch i = 0.6085rff (r = 2xlO-Vff ~ 3200(p,/102H2 cm-3)-i/2yr), 
Mcore ^ 0.4697cf/(47rG')=^/Vy^ ~ 0.11(c,/190ms-i)3(p,/102H2 cm-3)-V2Mo (model AHl) 
and Mcore ^ 0.41040^/(477(7)=^/^^^ - 0.09(c,/190ms-i)=^(p,/102H2 cm-3)-V2MQ (model 
AH2). From these results, it is shown that the core mass increases with time due to the 
continuous accretion and the mass does not depend on the exact equation of state in the 
core. This is understood as follows: the core mass is determined by the accretion rate of the 
isothermal gas which is independent from the polytropic F in the core. 

The gravity by the core has an effect on the outer inflow and outflow. Since the effect 
depends only on its mass, the difference of the polytropic index of the core does not play an 
important role for the inflow and outflow. Therefore, we will study this model AH to see the 
long time evolution of the outflow. 

In Figure 4c, the snapshot at t = 0.6105rff (r = 4 x 10~^Tff) is plotted for model AHl. 
Comparing this with Figure lb for model A (both have the same resolution but for different 
epochs), it is shown that the shock front which separates the inflow and the outflow passed 
the slow-mode MHD shock near z ^ O.Olif and has just reached the outer fast-mode shock 
front near z ~ 0.02H. The evolution of model AH2 is essentially the same. The maximum 
speed of the outflow reaches ~ 8cs ~ 2kms~^(cs/0.2kms~^). This maximum speed seems 
smaller than that observed in the molecular outflow. Since the mass accumulated in the 
core is only equal to ~ 0.1(c5/190ms~"^)^(ps/100H2 cm'~^)~^/^MQ, the outflow speed seems 
to be much faster than this value, when the mass has grown to typical T Tauri stars. Figure 
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3b indicates the outflow is accelerated near the core and the opening angle of the outflow- 
in this region is wide. However, departing from the acceleration region the flow changes its 
direction toward the z-axis. Figure 4c show that the opening angle decreases as the outflow 
proceeds. This indicates the flow is coUimated. 

3.4. Effect of the Initial Rotation Speed 

To see the effect of the initial rotation speed, we compare models AHl (Qq = Bt^^), BH 
(fio = 1%^), and CH (Vlo = 0.2rg^). These models have the same magnetic field strength, 
a. In Figure 5 panels (a-c), the structures at the final epoch of the isothermal run-away 
collapse phase are plotted for respective models. Models BH and CH indicate no prominent 
discontinuity in L6 (Figs. 5b and c), while model AHl has several shock fronts as described 
in §3.1. However, in LIO (not shown), there are discontinuities near 2; ~ 6 x 10~^if (model 
BH) and z ~ 5 x 10~^if (model CH) as weU as in model AHl (2; ~ 1 x 10~^i?). Comparing 
panels (b) and (c), distributions of the density and magnetic field lines are similar each other. 
This indicates that the evolution in the isothermal phase is slightly dependent on the initial 
angular momentum if ^ 1 x % ^■ 

Figures 5d-f show the structure at the age r = 4.5 x 10~^ts after the core formation. 
In all models the outflows are formed. However, the size of the region swept by the outflow 
is different for each model. With increasing Qq, more energetic outflow is driven. From the 
flow vectors, it is shown that model AHl (Fig.Gd) forms a bit more collimated outflow than 
models BH and CH (Fig.Ge and f). This seems to correspond to the differences in density 
distribution and magnetic field configuration. That is, in model AHl (also A and AH2) there 
is a relatively thick disk seen in L6 which is bounded by the shock fronts. This thick disk 
seems to confine the outfiowing gas in model AHl. While, in models BH and CH the disk 
is relatively thin, which seems to make the flow isotropic. Further, the opening angle of the 
magnetic fleld lines in models BH and CH is larger than that of model AHl. This causes 
the flow also open. 

Difference between models BH and CH comes from the fact that the epochs when the 
outflow begins are different. Since the initial angular momentum in model BH is flve-times 
larger than that of model CH, in model BH the outflow begins earlier than model CH. At 
T ~ 6 X 10~^Tff , however, even in model CH the top of the outflow reaches z ~ 0.02H and 
the structure looks very similar to model BH at r ~ 4.5 x 10~^Tff (Fig.5e). 

At the epoch when Figure 5d-f is taken (r = 4.5xlO~Vfr ~ 8xl0^yr[ps/10^H2 cm"^]"^/^), 
the mass in the adiabatic core reaches M^ore ^ 0.667cy{'inGf/'^pl^^ (model AHl), ~ 
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LSUcl/iAnGf/^pr (model BH), and ~ 1.469cf/(47rG')^/V^ (model CH). The instanta- 
neous rate of mass accretion to the adiabatic core for each model attains Mace = dMcorc/dt ~ 
110c3/(47rG') (model AHl), ~ 180c3/(47rG') (model BH), and ~ 220c3/(47rG) (model CH), 
respectively. Therefore, the core mass is approximately proportional to the mass accretion 
rate and the mass accretion rate is larger for models with smaller angular rotation speed flo. 

Accretion rate expected from the inside-out collapse model (Shu 1977) is equal to 
0.975cl/G = 12.25c^/{AnG). Therefore, the accretion rates calculated here are 9 — 18 times 
larger than that expected by the inside-out collapse model. 

Consider the reason why the mass of the core decreases with increasing Qq. Since the 
gas is supplied to the core mainly through the disk, we will consider the mass inflow/outflow 
transported in the disk. At that time, the gas disk can be divided into three regions. 
Outermost region is occupied with isothermal gas and the gas is contracting or inflowing 
(pseudo-disk). Therefore the inflow mass rate (Min) outermost > and the outflow mass rate 
(Mout) outermost = in this outermost region. Inside of this region, outflow is generated, 
although a large part of the gas is still inflowing. Therefore in this middle region, the inflow 
rate is smaller than that of the outermost region, (Min)middie ^ (Mn)outermost - (Mn)middie + 
(Mout)middie, ^ud the excess mass is transported to the outflow, (Mout)middie > 0. Innermost 
is the adiabatic core. Since the mass accretion rate to the core is equal to the net mass 
inflow rate from the middle region. Mace — ( Min) middle- 

Mass inflow driven by the self-gravity is more important in a model with small in 
which the self-gravity is ineffectively counterbalanced with the centrifugal force. Therefore 
(Min)outermost becomes larger for slow rotator. This is the first effect of the rotation. 

Furthermore, the outfiow brings away appreciable amount of gas. As mentioned pre- 
viously, the outfiow is strongly generated in the fast rotator. Thus, the mass outfiow rate 
increases with increasing ^Iq as (Mout)middic ~ 80c^/(47rG) (model AHl), ~ 20c^/(47rG) 
(model BH), and ~ 10c^/(47rG) (model CH). As a result, increasing Qq, the portion of out- 
flow gas to the inflow gas (Mout) middle /(Min) kiddie becomes large as ~ 40% for model AHl, 
~ 10% for model BH, and <^ 5% for model CH. These two effects works cooperatively to 
reduce the mass accretion rate Mace — (Min) kiddie to the core for the cloud with a large Qq. 

The maximum outflow speed realized in respective flgures are equal to Vmax — 9.3cs 
(model AHl), GAcg (model BH), and 3.2Cs (model CH). Since the outflow is accelerated by 
the toroidal magnetic fields which are generated by the rotation motion, this V^max increases 
with increasing Qq. 

As shown in Tomisaka (2000), since the excess angular momentum of the inflowing gas 
is effectively removed by the outflow, the total angular momentum of the first core which 
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is defined as a gas with p > pA is equal to Jcore — 9.2 x 10 ^CgH contained in the mass 
of M = 0.67cl/{AnGf/^py^ (model AHl), ~ 2.5 x IQ-^CsH in M = l.Slcl/iAnGf/^pl^^ 
(model BH), and ~ 7.4 x IQ-^CsH in M = 1.47^/ {AnGf/^pl^^ (model CH). The total 
angular momenta contained in the first cores are only 1.1%, 0.07%, and 0.9% of the initial 
ones contained in respective mass M. 

3.5. Effect of the Magnetic Field Strength 

To see the effect of the magnetic field strength, we calculated models NH, DH and EH 
in which we took a — 0, 0.1, and 0.01. 

3.5.1. Model without Magnetic Fields 

Model NH has no magnetic fields. In Figure 6a, a snapshot ait — 0.6977Tff captured by 
L8 is shown for model NH. At this stage, whole the cloud is in isothermal regime and the disk 
experiences the run-away collapse even if the centrifugal force may work to support the cloud. 
This confirms the earlier studies by Norman, Wilson, & Barton (1980) and Narita, Hayashi, & 
Miyama (1984). The physical reason why the centrifugal force docs not stop the contraction 
in the isothermal run- away collapse phase is explained in Hayashi (1980) as follows: Due 
to the the centrifugal force, the mass contained in the Jeans scale (~ Cg/ {Gpc)^^'^) from 
the center is decreasing throughout the collapse. In this sense the centrifugal force does 
work! Only small part of the mass that resides near the center becomes high-density. But 
contraction itself continues and the central density rises greatly in a finite time-scale, if the 
isothermal equation of state is valid. 

Similar to the previous magnetized models, a small adiabatic core is formed first. Since 
there is no magnetic fields, magnetic braking does not work, however, in this model. There- 
fore, gas that accretes on the core has relatively large angular momentum in contrast to the 
magnetized model. As a result, a ring forms by the gas which accreted on the adiabatic 
core. Since the specific angular momentum, j = rv^f,, increases further with time, radius of 
the ring grows. Another snapshot in panel (b) at t = 0.7011rfj (r ~ 3.4 x lO^^rg) shows the 
ring clearly. The ring seems unstable for non-axisymmetric perturbation. This may form 
a spiral structure similar to that found by Klein et al. (1999). However this is beyond the 
scope of this paper. Therefore, it is concluded that a rotating but non-magnetic cloud leads 
to a rotating ring in the adiabatic accretion stage. 
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3.5.2. Models with Magnetic Fields 

To see the effect of the magnetic field strength, in Figure 7 we compare models BH 
{a = 1), DH(a = 0.1), and EH {a = 0.01). All models have the same initial rotation speed 
Qq = 1, and polytropic gamma F = 2. In panels (a)-(c), the structures at the epoch when 
an adiabatic core begins to form are plotted. Comparing these panels, it is shown that in 
model B (a = 1) a flare-up disk is formed whose isodensity lines are departing from the disk 
mid-plane leaving from the center. Decreasing the initial magnetic field strength, the shape 
of dense part of the disk becomes rounder. Similar effect is already reported for non-rotating 
magnetized cloud collapse (Tomisaka 1996), that is, deceasing a the shape of the isothermal 
contracting disk becomes rounder and finally reaches a sphere for a = 0. 

In panel (d) we plotted a snapshot for model BH at t = 0.7264Tff (r = 4.46 x 10~^Tff) 
captured in L7. This is the same as illustrated in Figure 5c but captured in L7 which has 
twice as fine as Figure 5e. Figure 7d which shows the structure near the root of outflow 
indicates that it is very similar to that of model AHl. For example, the outflow leaves 
from the disk in the direction almost parallel to the disk but it changes its direction to the 
2;-direction. Time t = 0.7264rfj is equivalent to r = 4.46 x lO^^Tff. It is concluded that in 
a time-scale of r ~ 4.5 x lO^'^rgf, the flow pattern is completely changed from the run-away 
collapse to the outflow plus continuous inflow in the disk. The outflow gas flows along the 
surface whose shape resembles a capital letter U. The flow departs from the disk with a wide 
opening angle but it change its direction parallel to the 2;-axis. 

In panel (e) , we plotted the structure expected from a model with weaker magnetic fleld 
(model DH: a — 0.1 and flo — 1). The snapshot corresponds to the epoch oi t — 0.7836Tff. 
This corresponds to r = 4.59 x lO^^Tff which is similar to the time-scale between panels 
(a) and (d). In contrast to the previous model BH, the outflow gas forms a sphere and the 
magnetic field lines are folded inside this sphere. The magnetic field lines are folded by the 
pinch or hoop stress of the toroidal magnetic field. The toroidal magnetic field component 
is confined in the region where the adjacent poloidal magnetic field lines are running in 
the opposite directions. For example, the regions around (r, z) ~ (0.004if, 0.003if) and 
(r, z) ~ (0.006i/, O.OOSiJ) in Figure 7e. 

In this model, the initial poloidal magnetic fields are weak compared to model BH. 
Therefore, rotational motion amplifies the toroidal fields and their strength surpasses easily 
that of the poloidal ones. Thus, the hoop stress by the toroidal field pinches efficiently the 
magnetic fleld lines. In the outflow acceleration region, the toroidal component is predomi- 
nant over the poloidal one. Magnetic fleld lines are pinched locally and folded. As a result, 
a spherical magnetic bubble is formed in this process, in which the toroidal magnetic fleld 
is predominant. 
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Toroidal component of the magnetic fields is continuously generated by the twist motion 
driven by the disk rotation. The disk angular momentum is transferred by this process. As 
a result, we do not see the ring of model NH which is supported by the centrifugal force. 

For the model with extremely weak field, we calculated model EH {a — 0.01 and Q — 
1). In panel (f), we plotted the snapshot ai t — 0.7830Tff (r = 4.53 x 10~^Tff). Density 
distribution and magnetic field configuration show that the flow in the magnetic bubble is 
more complicated or turbulent in this model. The shape of the bubble is more elongated than 
that formed in model DH. Distribution of toroidal field lines does not show any systematic 
pattern inside the bubble. Size of the bubble both in the z- and r-directions are smaller 
than that of models BH and DH. It is concluded that the size of the outflow region increases 
with increasing the magnetic field strength a. Comparing these three models, it is shown 
that there are at least two types of outfiows. That is, a laminar U-type fiow in which fast 
moving gas fiows along the surface whose shape resembles a capital letter U and a turbulent 
outfiow in which the magnetic fields and the velocity are randomly oriented. 

The masses accumulated in r ~ 4.5 x 10~ Tfj are equal to 1.3cl/{ATTGf/^pl^'^ for model 
BH, 0.87c^/(47rG')=^/Vy^ for model DH, and 0.650^/(4776')=^/^^^ for model EH. This shows 
that the mass accretion rate. Mace, is an increasing function of the initial magnetic field 
strength, a. This seems strange if we remember Mace is an decreasing function of the initial 
rotation speed, Qq, since both a and Qq have an effect to counterbalance against the self- 
gravity. This means that the mass inffow rate in the isothermal run-away collapse region, 
(Mill) outermost; iucrcascs with increasing a. This seems to come from a number of reasons, that 
is, the initial cylindrical cloud becomes massive with increasing a. Another reason is related 
to the characteristic wave speed in the magnetized medium. That is, the characteristic speed 
of the fast mode MHD wave is equal to (cf + Bq/AttpY^^ — Cs{l + a)^/^ in the case that 
the wave is propagating perpendicular to the magnetic field lines. This implies that the 
mass infiow rates are proportional to (1 -|- cc)^/^ (A: B: C = 2.83: 1.15: 1.02). This is not 
inconsistent with the actual values. 



4. Discussion 

4.1. Evolution to Form the Second Core 

In the previous section, we introduced a simple polytrope above the critical density 
p > Pa except for model A. However, this is an approximation to sec a long evolution. In 
reality, the composite polytrope (eq. 14) should be used. Here, further evolution later than 
that shown in Figures 1 and 2 is shown. In model R, we used the composite polytrope 
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to include the combined effects such as the dynamical compressional heating, the radiative 
cooling through the dust thermal emissions, and the energy loss associated with the H2 
dissociation. This corresponds to the continuation of model A but we consider here a more 
compact cloud as Ps = 10*^H2 cm^'^. It should be addressed since ps is assumed hundred- 
times larger than other models, the size scale H is 10-times smaller than other models with 
Ps = 10^H2cm-3 as ~ 1.05 x 10iW(c,/190ms-i)(p,/10^H2 cm-3)-V2. 

The evolution up to the first core formation is similar to model A, that is, after the 
equation of state becomes hard, p > Pa, collapse of the adiabatic core slows down and 
isothermal gas begins to accrete on to the core. In this model since ps is assumed equal to 
10^H2 cm^'^, the size scale H is 10-times smaller than other models with ps = 10^H2 cm^'^. 
In Figure 8a, we plotted the structure captured by L6 just before the first core formation 
t — 0.720lTff. Since the accretion continues, the mass of the adiabatic core increases with 
time, which leads to collapse the core quasi-statically. The structure in this stage (at t — 
0.7239Tff or T = 3.8 x lO^^Tff) is shown in panel (b). The radius of the first core shows in 
this panel is equal to r ~ 3.5 x W'^H ~ 2.45(c,/190ms-i)(p,/10^H2cm~3)-i/2 ^p^-g 

quasistatic contraction phase ends when the central density reaches pe- The structure at 
this stage {t = 0.724230rff [r = 4.132 x lO^^rg]) is plotted in panel (c). Comparing this 
with panel (b), both of which illustrate L12, it is clear that the adiabatic core is contracting. 
After that, the second collapse begins. Since the thermal energy is lost by the process 
of the dissociation of H2, the equation of state is assumed soft (7 ~ 1.1) again. In this 
phase, the flow is very similar to that realized in the isothermal run- away collapse phase 
(p < Pa)- The similarity comes from the fact that the equation of states of the both 
phases are soft. In panel (d), we plotted the structure captured in L16 at t = 0.724236Tff 
(r = 4.138 X lO^'^Tff), which represents the typical structure in the second collapse phase. 
This continues until the central density exceeds p > pc, beyond which another adiabatic 
core (the second core) is formed shown in panel (e). The size of the second core is equal to 
~ 1.3x lO^^i^ ~ 2i?Q(cs/190ms"^)(p5/10^H2cm~^)~^/^. This meets a similar situation when 
the flrst adiabatic core is formed, that is, a central part of the gas obeys a harder equation 
of state and forms a quasistatic core, while the outer part obeys a softer equation of state 
and continues to collapse. Analogy between the first and the second core leads a bipolar 
outfiow. In panel (f), it is shown a second outfiow is formed around the second core. This 
is a snapshot of L16 at t — 0.724237Tff (r = 4.140 x 10~^Tff), which resembles the structure 
seen in Figure 3b. 

This model strongly indicates us that there is another kind of outfiow accelerated around 
the second core (the second outfiow) besides that formed around the first core (the first 
outfiow). The outfiow speed, ~ 50cs, is much faster than that of the first outfiow, due to the 
fast thermal speed in the adiabatic core and the deep local gravitational potential. Since the 
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time span of the simulation shown here is restricted, we can not trace the evolution further. 
However, the simulation predicts that at least two different outflows are formed each of 
which is related to different types of adiabatic cores. Since the flow speed of the second 
outflow is much faster than that of the first outflow, the respective outflows correspond to 
the molecular bipolar outflow (the first outflow) and the fast neutral wind or optical jets (the 
second outflow). The radial size of the outflow is approximately equal to ~ 2 x lO^^H ~ 



2.1 X 10i2cm(c,/190ms-i)(p,/10^H2cm-3)-i/2 (=o.l4 AU or 30Rq). This indicates that 



optical jets are found inside molecular bipolar outflows. 

In model R, wc assumed the ideal magnctohydrodynamics. However, the ionization 
fraction of high-density gas is quite low and the electric conductivity decreases as collapse 
proceeds. Calculation of ionization equilibrium (Nakano & Umebayashi 1986) shows us that 
after nn ^ lO^^cm"^ charged grains are more abundant than ions, as a carrier of electric 
charges. Since the mass-to-charge ratio of grains is much larger than that of ions, the electric 
conductivity decrease greatly and the magnetic field decays mainly through Joule dissipation. 
Therefore, in the late phase of the first core (quasistatic contraction), the magnetic field in 
the core decreases its strength till the field configuration becomes force-free. After the 
dissociation of H2 begins, the temperature of the central region is enough high to achieve the 
thermal ionization of metals. Thus, in the second collapse phase, coupling of magnetic fields 
is recovered. Prom these, the ideal MHD is consistent in the first and the second collapse 
phases. However, the magnetic fiux density at the beginning of the second collapse phases 
seems to be much weaker than that obtained here. As shown in §3.5.2, when the poloidal 
magnetic field is weak, the fiow pattern is different from the case that the magnetic energy 
is comparable to the thermal one. Therefore, if the magnetic fiux is partly lost from the 
central part of the first core, a turbulent outfiow around the second core seems to be formed 
similar to Figure 7f. 



Here in this subsection, mass infiow/outfiow rate and linear momentum outfiow rate 
are seen more closely. The outfiow mass loss rate through the boundary of each grid level is 
measured as 



where n represents the unit vector outwardly normal to the surface. The integrand means 
that only the outwardly running mass fiux is summed up. Similarly, the mass infiow rate is 



4.2. Mass Inflow/Outflow Rate 




boundary of Ln 




(26) 
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calculated as 

Mi„(Ln) = f pmax[-v ■n,0]dS. (27) 

J boundary of Ln 

These rates are calculated for respective levels of the nested grid system. 

The time variations of mass inflow and outflow rates are illustrated for respective models 
in Figure 9. The time (horizontal axis) is measured from the epoch of the core formation. 
Time variations, especially longer time variations, in inflow rates are common for all panels. 
Just after the core formation, the inflowing mass flux decreases with the distance from the 
center. This means that the inner part of the isothermal collapse region is well expressed by 
the Larson (1969)-Penston (1969) self-similar solution (Ogino, Tomisaka, & Nakamura 2000) 
which leads to the mass inflow rate of A7c^/G; Departing from the center, the approximation 
of the Larson- Penston solution becomes worse and the inflow becomes lower. At r ~ 4 x 
10~"^Tff the inflow rates for various surfaces converge on ~ 20c^/G, although intense outflow 
reduces the inflow rate below this value. 

We will see each model more closely. Outflow rate rises in the deeper levels (Ln with 
larger n) flrst and propagates to lower levels (Ln with smaller n). This indicates the outflow 
region expands both in r- and z-directions. Comparing models AHl (panel a), BH (panel 
b), and CH (panel c), the effects of the initial rotation speed, fl, are apparent. 

(1) with increasing Qq, the outflow begins earlier. 

(2) with increasing Qq, the mass outflow rate increases; although the ratio Mout/Afin is equal 
to only 10% in model CH {Qq = 0.2), in model AHl {Qq = 5) it attains ~ 50%. 

We compare models BH (panel b), DH (panel d), and EH (panel e) to see the effects 
of the initial magnetic fleld strength, a. This shows that with decreasing a (from models 
BH to DH) the mass outflow rate increases. However, in model E, in which extremely weak 
poloidal magnetic flelds is assumed, the time variation in the inflow and outflow rates is 
rather chaotic which is lead by the turbulent flow pattern realized in model EH. Averaging 
the rates as < M >= Mdt/T, the rate of model E is ~ 20 - 40% smaller than that of 
model DH. Considering the disk after rotating at a flxed angle, the disk of model EH can 
generate weaker toroidal magnetic flelds than model DH, since model EH has only weak 
poloidal (source) magnetic flelds. Therefore, it is reasonable that the mass outflow rate in 
model EH, < Mout >, is smaller than that of model DH. 

How about the increase from model BH to DH? This is inconsistent with the above 
discussion. This increase seems related to the fact that the flow patterns of models BH and 
DH are different. That is, in model DH a magnetic bubble, in which magnetic flcld lines 
are folded and amplifled, is formed and the bubble expands. In contrast, in model BH, the 
gas is flowing outward along speciflc magnetic fleld lines and outflowing gas forms a wall of 
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a chimney whose shape looks hke a capital letter U. Mass outflow rate seems larger in the 
magnetic bubble type outflow rather than the chimney wall type. 

The linear momentum outflow rates are calculated by 



The time variation of momentum outflow rate is shown in Figure 10. 

Panels (a), (b), and (c) have a prominent common feature that the momentum outflow 
rates (dMV/dt) calculated at the boundaries of Ln {n < 7) are much larger than that of 
the boundaries of Ln {n > 7). This feature is common in models with a — 1 (models 
AHl, AH2, BH, and CH). However, models with weaker magnetic field does not show this 
feature. Recall the fact that models AHl, AH2, BH, and CH form a U-shaped outflow and 
in contrast models DH and EH form turbulent magnetic bubble type outflow. Difference 
in the momentum outflow rate seems to be related to the outflow pattern. In the laminar 
U-shaped outflow, the outflow is ejected with a rather wide opening angle and coUimated 
to the symmetric axis. Since equation (28) counts the momentum passing the upper and 
lower boundaries, the outflow rate increases after the outflow is coUimated and the gas flows 
parallel to the 2;-axis. This occurs in the larger scale than L7. 

Typical momentum outflow rate for molecular bipolar outflows observed by ^^CO(J = 
2 — 1) (Bontemps et al. 1996) are distributed between ~ lO~"^M0yr~-^ kms~^ for objects 
associated with active Class IR sources and ~ 2 x lO~^M0yr~^ kms~^ for objects associated 
with late Class 1 IR sources. This range corresponds to 320cf/G — 6.5cj/G, respectively, if 
we assume = 190m s~^. Figure 10 shows us that the momentum outflow rate in the range 
from ~ lOCg/G to ~ AOCg/G is expected. Thus, the momentum outflow rate for almost all 
the CO bipolar outflow sources associated with low-mass young stellar objects except for the 
active early class sources are explained by this model. 

Comparing models AHl, BH, and CH, it is shown that the momentum outflow rate 
increases with Qq- Comparing mass outflow and momentum outflow rates, models AHl, 
BH, and CH indicate that the momentum outflow rate is approximately proportional to the 
mass outflow rate, which means that the outflow speed is approximately equal irrespective of 
Qq. However, changing a is more complicated. Compare panels (b), (d), and (e). As shown 
in Figure 9, model DH {a = 0.1) shows larger mass outflow rate than model BH (a = 1). 
In contrast, as for the momentum outflow rate, the maximum rate of model BH is larger 
than that of model DH (Fig. 10). This represents the outflow speed of model BH {a — 1) is 
faster than model DH (a — 0.1). Panel (e) shows that the momentum outflow rate has a 
chaotic time variation as well as the mass outflow rate (Fig.9e), for model EH, which seems 
to correspond the turbulent outflow shown in Figure 7f. 



dMV 
dt 



out 




upper and lower boundaries of Ln 



pvz max[v • n, 0]dS. 



(28) 
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Summarizing the effect of changing a, 

(1) the mass outflow rate increases with increasing a ^ 0.1 but decreases after a ~ 0.1. 

(2) the outflow speed is an increasing function of a. 

(3) there are distinctly different two types of flow patterns: models with strong magnetic 
fields lead the U-shape outflow, while the flow becomes turbulent for models with weak 
magnetic fields. 



4.3. Where Is the Gas Ejected? 

Here, we intend to clear where the outfiow is ejected. In the run-away collapse phase, 
the rotational motion is relatively unimportant. After the core formation, rotation motion 
increases especially around the first core. This occurs near the centrifugal radius, which is 
defined for a gas with specific angular momentum j as 

R. - Ce^, (29) 

where M denotes the mass inside the radius Rc and Cc represents a numerical factor of the 
order of unity. Considering the self-similar solution for the run-away collapse phase (Saigo & 
Hanawa 1998), the specific angular momentum j is proportional to M. Further a uniform- 
density cylinder which is rotating with a uniform rotation speed u has a j distribution 
proportional to M. Assuming j = qGM/cg {q is a numerical factor), the centrifugal radius 
is written 

= c^- ■ (30) 

As shown in section 3, the angular momentum is redistributed in one magnetic flux tube. 
That is, gas near the disk surface obtains large amount of the specific angular momentum 
but that near the disk mid-plane loses the angular momentum. If the gas near the disk 
surface has the angular momentum of j+ > j, the condition that the effective potential at 
the centrifugal radius 

*(^) = -f-|f (^^) 

is larger than or equal to zero leads to a minimum angular momentum of 



J+ > (2Ce)V^j. 



(32) 
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5. Summciry 

We have explored the evolution of a magnetized interstellar cloud rotating around the 
symmetric axis. Following the change in the equation of state of the interstellar gas (Tohline 
1982), the cloud experiences several phases before going to a star, that is, the isothermal 
run-away collapse, the slowly contracting core composed of the molecular hydrogen (the 
first core), the second run-away collapse in the high-density gas where the dissociation of 
hydrogen molecules proceeds, and finally the second core which is made up of the atomic 
hydrogen. The magnetized cloud forms a pseudo-disk in these first and second run-away 
collapses, in which a supersonically contracting disk is formed and magnetic field lines are 
running perpendicularly to the disk. In the pseudo-disk, a number of fast- and slow-mode 
MHD shock pairs are formed which is extending parallelly to the disk. Just after the core 
is formed at the center, an accretion shock front appears which surrounds the core, through 
which the supersonic inflow motion is decelerated. While the flrst and second cores are 
slowly contacting, the outer pseudo disk continues to contract. Just outside the accretion 
shock front, the infall motion is accelerated and thus rotational motion becomes important 
from the conservation of angular momentum. By the effect of rotational motion, the toroidal 
magnetic fields and the poloidal currents are amplified, which bring a strong magnetic torque. 
The magnetic torque leads the angular momentum transfer from mid-plane to surface of the 
disk. This is the origin of the outflow found in star forming regions. Large scale bipolar 
molecular outflows are made outside of the flrst core, while optical jets and fast neutral 
winds are expected to be accelerated outside of the second core. Matter lost its excess 
angular momentum continues to contract further to form a star. 

This work was supported partially by the Grants-in-Aid (11640231, 10147105) from 
MEXT (the Ministry of Education, Culture, Sports, Science and Technology). Numerical 
calculations were carried out by Fujitsu VPP5000 at the Astronomical Data Analysis Center, 
the National Astronomical Observatory, an inter-university research institute of astronomy 
operated by MEXT, Japan. 
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Table Captions 



Table 1: Conversion from the non-dimensional to the physical quantities . 



Physical Quantities 


Conversion Factors 


Physical Values (// = 2.33) 


Velocity 


Cs 


190m s-i 


Density 


Ps 


lO^Hs cm-3 


Length 


c,/(47rG'p,)V2 ^ H 


0.341(c,/190ms-i)(p,/102H2 cm-3)-V2pc 


Time 


(47rG'p,)-i/2 = Tff 


1.75(p,/102H2 cm-3)-V2Myr 


Mass 




0.227(c,/190ms-i)3(p,/102H2 cm-'^)-^''^MQ 


Mass accretion rate 




1.29 X 10-^(c,/190ms-i)3Moyr-i 




cl/G- 


1.62 X 10-6(c5/190ms-^)3MQyr-i 


Momentum inflow/outflow rate 


ctlAiiG 


2.45 X lO-8(c,/19Oms-i)^M0yr-ikms"^ 




4IG^ 


3.08 X lO-^(c,/19Oms-i)4M0yr-^kms-^ 


Magnetic Field 


1/2 
Csps 


3.75(c,/190m s"^) (p./lO^Ha cm'^) V2^g 



"To meet the conventional normalization, we adopt c^/G in section 4. 
''To meet the conventional normalization, we adopt /G in section 4. 
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Table 2: Model Parameters. 



Model 


a 


Qo 


Pa 


Ps 


polytrope 










(H2 ciii-'^O 




A ... 


1 


5 


10^ 


102 


Realistic 


AHl . . . 


1 


5 


108 


102 


r = 2 


AH2 . . . 


1 


5 


10^ 


102 


r = 5/3 


B ... 


1 


1 


10^ 


102 


Realistic 


BH ... 


1 


1 


10^ 


102 


r = 2 


CH ... 


1 


0.2 


10^ 


102 


r = 2 


DH ... 


0.1 


1 


108 


102 


r = 2 


EH ... 


0.01 


1 


108 


102 


r = 2 


N ... 





5 


108 


102 


Realistic 


NH ... 





5 


108 


102 


r = 2 


R ... 


1 


1 


10^ 


10^ 


Realistic 
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(a) 
(b) 

(c) 

Fig. 1. — Evolution of model A with a = 1 and Qq = 5. Snapshots at the time of t = 0.6066rgf 
captured by different levels are shown: LI (a), L5 (b), and LIO (c). The actual size of the 
frames of L5 (b) and LIO (c) are, respectively, 1/16 and 1/512 smaller than that of LI (a). 
Magnetic field lines (almost vertical) and isodensity contours are presented as well as the 
velocity vectors by arrows. Near the right-bottom corner, the logarithm of the maximum and 
the minimum of the densities are numerically shown. The maximum speed is also shown 
and the velocity vector corresponding its value is illustrated at the lower-left corner by a 
horizontal arrow. 



Fig. 2. — Cross-cut views along the equatorial plane (a) and the 2;-axis (b). The figures 
show formation of the first core. In panel (a), log p{r,z — 0) (sohd fines), log Bz{r,z — 0) 
(short dashed lines), v^{r,z — 0) (long dashed lines), and —Vr{r,z — 0) (dotted lines) are 
plotted, while in panel (b), logp(r = 0, z) (solid lines), logBz{r — 0, z) (short dashed lines), 
and —Vz{r — 0, z) (dotted lines) are plotted. 

(a) 

(b) 
(c) 

Fig. 3. — The same as Fig.l but for the snapshot at t = 0.6069rff. Panel (a) shows the 
structure captured by LIO which is the same as Fig.lc. At this stage, gas begins to outflow 
from the disk. Outflow sweeps the sphere r ^ 1.2 x IQ^^H . Panel (b) corresponds to L12, 
which has 4 times flner spatial resolution than panel (a). Panel (c) shows toroidal-to-poloidal 
magnetic field strength ratio. Contour lines extending in the direction ±45 deg from the disk 
are for the toroidal-to-poloidal ratios. Poloidal magnetic field fines are the same as panel (b). 
It is shown that the gas outfiows along the region where toroidal magnetic field is important. 

(a) 
(b) 
(c) 

Fig. 4. — The same as Fig.l but for models AHl and AH2. Snapshots at the same epoch 
of Fig.3a, t = 0.6069rff, are shown for models AHl (F = 2) and AH2 (F = 5/3) in panels 
(a) and (b), respectively. Although the structure of the core is different, the outfiow is very 
similar with each other. In panel (c), snapshot at t = 0.6105rff (r = 3.94 x 10~^rfj from 
the core formation) is plotted for model AHl. Be careful that the linear size of this panel 
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is 32 times larger than panels (a) and (b). Comparing with Fig. lb (the same resolution), it 
is shown that the shock front passed the slow-mode MHD shock and has just reached the 
outer fast-mode shock front. 

(a) (b) 
(c) (d) 
(e) (f) 
Fig. 5. — Comparison of models with the same magnetic field strength, a, but different 
rotation speed, Qq. Panels (a) and (d) represents the structure of model AHl (Oq = ^t^^) 
captured by L6. Panels (b) and (e) are for model BH. Model B corresponds to slower 
rotator Qq = Ir^^. Panels (b) and (e) represents the structure captured by L6 at the ages 
of t = 0.7219rfr and t = 0.7264rff (r = 4.5 x W~^ts), respectively. Panels (c) and (f) are 
for model CH (Oq = 0.2rff) and show the snapshots ai t — 0.7262rff and t — 0.7307rff 
(r = 4.5 X 10~^Tff), respectively. 

(a) (b) 
Fig. 6. — The same as Fig.l but for model NH. This model is for a non-magnetized cloud. 
In panel (a), a snapshot at t = 0.6977rff captured by L8 is shown. At this stage, whole the 
cloud is in isothermal regime. Another snapshot at t — 0.701lTff is shown in panel (b). Gas 
that accretes afterwards forms a ring which is supported essentially by the centrifugal force. 

(a) (b) 
(c) (d) 
(e) (f) 
Fig. 7. — Comparison of models with the same rotation speed D,o — 1 but different magnetic 
field strength a. Panels (a), (b), and (c) represent the structure when the adiabatic core 
begins to form, while panels (d), (e), and (f) represents that after a protostar is formed. In 
panel (a) we plotted a snapshot for model BH captured by L7 at t — 0.7219rff, which is the 
same snapshot shown in Fig. 5b but for different level. Panel (d) represents the structure of 
the same model but for the protostellar phase, that is, t = 0.7264rff (r = 4.46 x lO'^Tfj). 
This corresponds to Figure 5e. Panels (b) and (e) illustrate snapshots at t = 0.7790rfj and 
that at t = 0.7836rff (r = 4.59 x IQ-Vg), respectively, for model DH (a = 0.1). Panels (c) 
and (f) are for a model with extremely weak magnetic fields (model EH; a — 0.01). The 
snapshots at the epoch t — 0.7784Tff and t — 0.7830Tff (r = 4.53 x 10~^Tfr) are illustrated in 
panels (c) and (f), respectively. 

(a) (b) 
(c) (d) 
(e) (f) 
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Fig. 8. — The evolution of model R. In this model since ps is assumed equal to 10^H2 cm~^, 
the size scale H is 10-times smaller than other models with ps = 10^H2cm~^. In panel (a), 
we plotted the structure captured by L6 just before the first core formation t = 0.7201rff . At 
t = 0.7239rfj (r = 3.8 x 10~^rff), the first core gradually contracts by the effect of continuous 
mass accretion (panel b). Finally, t = 0.724230rfr (r = 4.132 x 10~^rfr), the central density 
reaches pe (panel c). After that, the second collapse begins. In this phase, flow is very 
similar to that realized in the isothermal run-away collapse phase. In panel (d), we plotted 
the structure captured in L16 at i = 0.724236rff (r = 4.138 x 10~^Tfr). After the central 
density exceeds pc, another adiabatic core (the second core) is formed (panel e). By a similar 
mechanism to form a bipolar outflow, a second outflow is formed around the second core 
(panel f). This is a snapshot of L16 at t = 0.724237Tfr (r = 4.140 x IO'^ts). 



Fig. 9. — Mass inflow/outflow rates measured at the boundaries of nested grid systems 
(eqs.27 and 26). Panel (a)-(e) corresponds to models AHl, BH, CH, DH, and EH. The 
horizontal axis represents the time after the core formation. Levels at which the mass 
inflow/outflow rates are measured are shown near the respective lines. We plotted Lll in a 
solid line, LIO in a dotted line, L9 in a short-dashed line, L8 in a long-dashed line, L7 in a 
short dash-dotted line, and L6 in a long dash-dotted line. In panel (e) to specify which line 
represents the outflow, we added brackets to the outflow rates. Inflow rates are larger than 
outflow rates for respective levels. 



Fig. 10. — Linear momentum outflow rates measured at the upper and lower boundaries of 
nested grid systems (eqs.28). Panel (a)-(e) corresponds to models AHl, BH, CH, DH, and 
EH. The horizontal axis represents the time after the core formation. Levels at which the 

linear momentum outflow rates arc measured are shown near the respective lines. We plotted 
Lll in a solid line, LIO in a dotted line, L9 in a short-dashed line, L8 in a long-dashed line, 
L7 in a short dash-dotted line, and L6 in a long dash-dotted line. 
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